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p^ ' Abstract 

_^ In this paper, a parametric level set method for reconstruction of obsta- 

cles in general inverse problems is considered. General evolution equations 
for the reconstruction of unknown obstacles are derived in terms of the un- 
derlying level set parameters. We show that using the appropriate form 
of parameterizing the level set function results a significantly lower dimen- 
sional problem, which bypasses many difficulties with traditional level set 
\^ methods, such as regularization, re-initialization and use of signed distance 

r^^ function. Moreover, we show that from a computational point of view, low 

^O order representation of the problem paves the path for easier use of New- 

ton and quasi- Newton methods. Specifically for the purposes of this paper, 
we parameterize the level set function in terms of adaptive compactly sup- 
ported radial basis functions, which used in the proposed manner provides 
(^ flexibility in presenting a larger class of shapes with fewer terms. Also they 

^^ provide a "narrow-banding" advantage which can further reduce the number 

^ of active unknowns at each step of the evolution. The performance of the 

k> proposed approach is examined in three examples of inverse problems, i.e., 

5—1 electrical resistance tomography. X-ray computed tomography and diffuse 

optical tomography. 
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1 Introduction 



Inverse problems arise in many applications of science and engineering including 
e.g., geophysics 67,85 , medical imaging [3,48,80 , nondestructive evaluation 47 
49 and hydrology 



13 71,82 



In all cases the fundamental problem is usually 
extracting the information concerning the internal structure of a medium based on 
indirect observations collected at the periphery, where the data and the unknown 
are linked via a physical model of the sensing modality. A fundamental challenge 
associated with many inverse problems is ill-posedness (or ill-conditioning in the 
discrete case), meaning that the solution to the problem is highly sensitive to noise 
in the data or effects not captured by the physical model of the sensor. This 
difficulty may arise due to the underlying physics of the sensing system which in 
many cases (e.g., electrical impedance tomography [ITj, diffuse optical tomography 
[3], inverse scattering [24], etc) causes the data to be inherently insensitive to 
fine scale variations in the medium. This phenomenon makes such characteristics 
difficult, if not impossible to recover stably 74 . Another important factor causing 



the ill-posedness is limitations in the distribution of the sensors yielding sparse data 
sets that again do not support the recovery of fine scale information uniformly in 
the region of interest 52 . Many inverse problems of practical interest in fact suffer 
from both of these problems. From a practical point of view, left untreated, ill- 
posedness yields reconstructions contaminated by high frequency, large amplitude 
artifacts. 

Coping with the ill-posedness is usually addressed through the use of regulariza- 
Based on prior information about the unknowns, the regularization 
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tion 

schemes add constraints to the inverse problem to stabilize the reconstructions. 
When the inverse problem is cast in a variational framework, these regularization 
methods often take the form of additive terms within the associated cost function 
and are interpreted as penalties associated with undesirable characteristics of the 
reconstruction. They may appear in various forms such as imposing boundedness 
on the values of the unknown quantity (e.g., Tikhonov or minimum norm regular- 
izations 30,74]) or penalizing the complexity by adding smoothness terms (e.g.. 



the total variation regularization [I]). These regularization schemes are employed 
in cases where one seeks to use the data to determine values for a collection of un- 
knowns associated with a dense discretization of the medium (e.g, pixels, voxels, 
coefficients in a finite element representation of the unknown [63]). 

For many problems, the fundamental objective of the process is the identi- 
fication and characterization of regions of interest in a medium (tumors in the 



body ll2j, contaminant pools in the earth 34 , cracks in a material sample pOl, etc). 



For such problems, an alternative to forming an image and then post-processing 
to identify the region is to use the data to directly estimate the geometry of the 
regions as well as the contrast of the unknown in these regions. Problems tack- 
led in this way are known as the inverse obstacle or shape-based problems. For 



earlier works in this area see 18,40,43 and particularly the more theoretical ef- 
forts by Kirsch I39j and Kress et al. |41| . Such processes usually involve a rather 
simple parametrization of the shape and perform the inversion based on using the 
domain derivatives mapping the scattering obstacle to the observation. Relative 
to pixel-type approaches, these geometric formulations tend to result in better 
posed problems due to a more appropriate obstacle representation. Moreover, in 
such problems the regularization can be either performed implicitly through the 
parametrization or expressed in terms of geometric constraints on the shape 
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However, this class of shape representation is not topologically flexible and the 



number of components for the shape should be a priori known 64 . They also 
create difficulties in encountering holes and high curvature regions such as the 
corners. These difficulties have lead over the past decade or so to the develop- 
ment of shape-based inverse methods employing level set-type representation of 
the unknowns. 



The concept of level sets was first introduced by Osher and Sethian in 57 
This method was initially designed for tracking the motion of a front whose speed 
depends on the local curvature. The application of the level set approach to inverse 



problems involving obstacles was discussed by Santosa in 64 . One of the most 



attractive features of the level set method is its ability to track the motion through 
topological changes. More specifically, an a priori assumption about the connect- 
edness of the shapes of interest was no longer required. Following Santosa's work. 



Litman et al. in 46 explored the reconstruction of the cross-sectional contour of 



a cylindrical target. Two key distinguishing points about this work were in the 
way that the authors dealt with the deformation of the contour, and in their use 
of the level set method to represent the contour. The shape deformation method 
implemented in this work was enabled by a velocity term, and lead to a closed-form 
derivative of the cost functional with respect to a perturbation of the geometry. 
They defined shape optimization as finding a geometry that minimizes the error in 



the data fit. Later Dorn et al. in 26 introduced a two-step shape reconstruction 
method that was based on the adjoint field and level set methods. The first step 
of this algorithm was used as an initialization step for the second step and was 
mainly designed to deal with the non-linearities present in the model. The second 
step of this algorithm used a combination of the level set and the adjoint field 
methods. Although inspired by the works of [46||57|[64] , the level set method used 



by Dorn et al. was not based on a Hamilton- Jacobi type equation, instead, an opti- 
mization approach was employed, and an inversion routine was applied for solving 
the optimization. The level set ideas in inverse problems were further developed 
to tackle more advanced problems such as having shapes with textures or multiple 
possible phases 15,25 . Moreover, regarding the evolution of the level set function 



where usually gradient descent methods are the main minimization schemes ap- 
plied, some authors such as Burger 10 and Soleimani 68 proposed using second 
order convergent methods such as the Newton and quasi-Newton methods. 

Although level set methods provide large degrees of flexibility in shape repre- 
sentation, there are numerical concerns associated with these methods. Gradient 
descent methods used in these problems usually require a long evolution process. 
Although this problem may be overcome using second order methods, the per- 
formance of these methods for large problems such as 3D shape reconstructions 
remains limited and usually gradient descent type methods are the only option 
for such problems. Moreover re-initialization of the level set function to keep it 
well behaved and velocity field extensions to globally update the level set function 
through the evolution are usually inevitable and add extra computational costs 
and complexity to the problem 56 . More detailed reviews of level set methods in 



inverse problems can be found in 11,24 



In all traditional level set methods already stated, the unknown level set func- 
tion belongs to an infinite dimensional function space. From an implementation 
perspective, this requires the discretization of the level set function onto a dense 
collection of nodes. An alternative to this approach is to consider a finite dimen- 
sional function space or a parametric form for the level set function such as the 
space spanned by a set of basis functions. Initially, Kilmer et al. in 38 proposed 
using a polynomial basis for this purpose in diffuse optical tomography applica- 
tions. In this approach, the level set function is expressed in terms of a fixed order 
polynomial and evolved through updating the polynomial coefficients at every it- 
eration. Parametrization of the level set function later motivated some authors in 
field of mechanics to use it in applications such as structural topology optimiza- 
One of the main contributions in this regard is the work by Wang et 



tion 
al. in 
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Here the level set function is spanned by multiquadric radial basis func- 
tions as a typical basis set in scattered data fitting and interpolation applications. 
Authors showed that through this representation, the Hamilton- Jacobi partial dif- 
ferential equation changes into a system of ordinary differential equations and the 
updates for the expansion coefficients may be obtained by numerically solving an 
interpolation problem at every iteration. 

More recently, the idea of parametric representation of the level set function has 



been considered for image processing applications |5j[29j. Gelas et al. in 29 used 
a similar approach as the one by Wang et al. for image segmentation. As the basis 
set they used compactly supported radial basis functions, which not only reduce the 
dimensionality of the problem due to the parametric representation, but also reduce 
the computation cost through the sparsity that this class of functions provide. 
As advantages of the method they showed that appropriate constraints on the 
underlying parameters can avoid implementation of the usual re-initialization. Also 
the smoothness of the solution is guaranteed through the intrinsic smoothness of the 
underlying basis functions, and in practice no further geometric constraints need 
to be added to the cost function. As an alternative to this approach, Bernard et al. 
in [5] parameterized the level set function in terms of B-splines. One of the main 
advantages of their method was representing the the cost function minimization 
directly in terms of the B-spline coefficients and avoiding the evolution through 
the Hamilton- Jacobi equation. 

In this paper the general parametric level set approach for inverse problems is 
considered. For an arbitrary parametrization of the level set function, the evolution 
equation is derived in terms of the underlying parameters. Since one of the main 
advantages of this approach is low order representation of the level set function, 
in practice the number of unknown parameters in the problem is much less than 
the number of pixels (voxels) in a traditional Hamilton- Jacobi type of level set 
method, therefore we concentrate on faster converging optimization methods such 
as the Newton or quasi-Newton type method. To represent the parametric level 



set function, as in 29 , 79 we have proposed using radial basis functions. However 



unlike the previous works employing radial basis functions in a level set frame- 
work, in addition to the weighting coefficients, the representation is adaptive with 
respect to the centers and the scaling of the underlying radial basis functions. This 
technique basically prevents using a large number of basis terms in case that no 
prior information about the shape is available. To fully benefit from this adaptiv- 
ity, we further narrow our choice by considering compactly supported radial basis 
functions. Apparently, this choice would result sparsity in the resulting matrices. 
However, we will discuss a behavior of these functions which can be exploited to 
further reduce the number of underlying basis terms and provide the potential 
to reconstruct rather high curvature regions. The flexibility and performance of 
proposed methods will be examined through illustrative examples. 

The paper is structured as follows. In Section |2] we review shape-based inverse 
problems in a general variational framework. Section |3] is concerned with obtaining 
the flrst and second order sensitivities of the cost function with respect to the 
functions deflning the shape. Based on the details provided, a brief revision of the 



relevant traditional level set methods is provided paving the path to use parametric 
level set methods. In Section |4] the general parametric level set method will be 
introduced and evolution equations corresponding to the underlying parameters 
are derived. In Section [5] an adaptive parametric level set based on the radial basis 
functions will be proposed and the approach will be narrowed down to compactly 
supported class of functions, due to their interesting properties. Section [6] will 
examine the method through some examples in the context of electrical resistance 
tomography, X-ray tomography and diffuse optical tomography. Finally in Section 
[7] some concluding remarks and suggestions will be provided. 

2 Problem Formulation 

2.1 Forward Modelling 

The approach to modelling and determination of shape we consider in this paper 
is rather general with the capability of being applied to a range of inverse prob- 
lems. In Section |6] we specifically consider three applications: electrical resistance 
tomography, limited view X ray tomography, and diffuse optical tomography. The 
details of the specific problems will be provided in the same section within the 
context of the examples themselves. Up until that point, we have chosen to keep 
the discussion general. 

Consider f2 to be a compact domain in M", n > 2, with Lipschitz boundary 
dQ. Further assume for x € i7, a space dependent property p(x) e Sp, where 
Sp is a Hilbert space. A physical model Ai acts on a property of the medium 
(e.g., electrical conductivity, mass density or optical absorption), p, to generate an 
observation (or measurement) vector u, 

u = M{p), (1) 

where u itself belongs to some Hilbert space S^. In most applications m is a vector 
in C', the space of k dimensional complex numbers where k represents the number 
of measurements and accordingly a canonical inner product is used. 

As a convention throughout this paper, to keep the generality of notation, the 
inner products and norms corresponding to any Hilbert space H, are subindexed 
with the notation of the space itself, e.g. (., .)h. 



2.2 Inverse Problem 

The goal of an inverse problem is the recovery of information about the property p 
based on the data u. Here we consider a variational approach where the estimate 
of p is generated as the solution to an optimization problem. The functional 
underlying the problem is usually comprised of two terms. The first term demands 
that the estimate of p be consistent with the data in a mathematically precise 
sense. As is well known however, many interesting inverse problems are quite ill- 
posed. This means when the data consistency is our only concern, the resulting 
estimate of p could be quite far from the truth, corrupted by artifacts such as high 
frequency and large amplitude oscillations. Hence, an additional term (or terms) 
are required in the formulation of the variation problem which capture our prior 
knowledge concerning the expected behavior of the p in Q. Such terms serve to 
stabilize (or regularize) the inverse problem. Defining a residual operator as 

n{p) = M{p)-u, (2) 

the inverse problem is formulated in the following manner 

mm ^{p) = ^-\\n{p)\\l+^{p), (3) 

where ^ is the regularization functional. Appropriate choice of Jtf is usually based 
on properties of the problem. In the typical case where the unknown property p is 
represented as a dense collection of pixels (or voxels) in Q, the regularization penal- 
ties are used to enforce smoothness and boundedness of p in Q. These are usually 
considered in the framework of Tikhonov and total variation regularizations [1,74 



An alternative approach that has been of great interest in recent years is based 
on geometric parameterizations of the unknown. Here the regularization penalties 
are either embedded in the nature of the unknown or expressed as geometric con- 
straints on the unknown. No matter which approach is used, usually in defining 
^ different spaces and their corresponding norms may be used. For a more de- 



tailed review of such methods an interested reader is referred to 24 , 77 . Clearly 
when the number of parameters involved in the problem is sufficiently small, an 
underdetermined inverse problem can be made overdetermined, and in this sense 
the problem becomes better posed. Since the parametrization idea we will put 
forth in this paper is empirically found to be well-posed enough that no necessary 
regularization terms need to be added to the cost function, ^ will be neglected in 
our future discussions of ^. 



2.3 A Shape-Based Approach For the Unknown Parameter 



For a large class of "shape-based" inverse problems I24l^2lpllm4] , it is natural to 
view p(x) as being comprised of two classes, i.e., foreground and background. The 
problem then amounts to determination of the boundary separating these classes as 
well as characteristics of the property values in each class. In the simplest case, p(x) 
is piecewise constant while in more sophisticated cases p(x) may be random and 
characterized by different probabilistic models in the two regions. In this paper, 
we assume that over each region p{x.) is at least differentiable. The property 
of interest in this case is usually formulated through the use of a characteristic 
function. Given a closed domain D '^ Q with corresponding boundary dD, the 
characteristic function xd is defined as 



Accordingly, the unknown property p(x) can be defined over the entire domain fl 
as 

p(x) =pi(x)xD(x)+Po(x)(l-XD(x)). (5) 

Unlike p(x) which is clearly not differentiable along dD, Pi(x) indicating the prop- 
erty values inside D and Po(x) denoting the values outside, are assumed to be 
smooth functions, and as mentioned earlier, at least belonging to C^{Q). 

In a shape-based approach, finding dD is a major objective. As ([s]) and Q 
show, p(x) is implicitly related to D and to find dD, a more explicit way of relating 
them should be considered. In this regard the idea of using a level set function 
proves to be especially useful ||64j. Here dD is represented as some level set of a 
Lipschitz continuous function : f2 ^ M. When the zero level set is considered, 
0(x) is related to D and dD via 

(/)(x) > Vx e D 

0(x) = y^edD (6) 

0(x) < Vx € fi \ D. 

Making the use of a Heaviside function, defined as H{.) = ^(l + sign(.)), the 
function p(x) can be represented as 

p(x) = p,(x)if (0(x)) + p„(x)(l - i/(0(x))). (7) 

This equation in fact maps the space of unknown regions D into the space of 
unknown smooth functions (h. 



3 Inversion as a Cost Function Minimization 

In this section we develop the mathematical details of the minimization problem ^ 
when a shape-based approach as ([T]) is considered. Most current methods use the 
first and second order sensitivities of the cost function with respect to the unknown 



parameter to perform the minimization 11 . Based on the general details provided 



we will briefly revisit the traditional level set approaches relevant to this paper in 



Section 3.2, since understanding the details better justifies the use of parametric 



level set methods. 

3.1 Cost Function Variations Due to the Unknowns 

We begin by assuming that the first and second order Frechet derivatives of Tl{p) 
exist and denote them as 'JZ'{p)[ ■ ] and 7l"{p)[., .]. The first order Frechet deriva- 
tive of a function (if it exists) is a bounded and linear operator. The second order 
Frechet derivative is also bounded but bilinear, which means the operator acts on 
two arguments and is linear with respect to each El. 

For an arbitrary variation 6p e Sp and the real scalar e, using the generalized 
Taylor expansion we have 

n{p + eSp) = n{p) + en'{p)[5p] + —n"{p)[8p, Sp] + Oie"^). (8) 

Rewriting ^{p) as 

^ip) = l{np),np)h.. (9) 



and recalling the fact that {ui,U2)s^ = (^2, ^i)§„ for ^i, ^2 ^ S„ and overline denoting 
complex conjugate, the variations of the cost function with respect to the variations 
of p can be derived as 

^{p + e6p) = ^{p) + e^'{p)[6p\ + —^"{p)[5p, 6p] + 0{e^), (10) 

where for pi,p2 ^ Sp 

^\p)[pi]=Me{n\p)[p,mp))^^ (11) 

and 

=^"b)bi,P2] =■^e(7^'(p)bl],7^'(p)b2])§. +i^e(7^"(p)bl,p2],7^(p))s„. (12) 



The notation ^e indicates the real part of the corresponding quantity. Denoting 
TZ'ipyi ■ ] ^s the adjoint operator between §„ and Sp as 

{u,n'(p)[p])s^ = {n'ipyiuiph,, vn e §„, vp € Sp, (13) 



(11) can be written as 

^'{p)[pi] = Me{n'{py[nip)ip,h,. (u) 



Equations (11), (14) and (Il2j) are in fact the first and second order Frechet deriva- 
tives of ^ with respect to p. In a more general context (and indeed one which we 
shall use in Section 111), p itself can be the map over some variable v from another 
Hilbert space S„ into Sp, i.e., p(v) : §„ -^ Sp. Assuming the existence of the first 
and second order Frechet derivatives of p with respect to v, denoted as p'{v)[ . ] 
and p"{v)[., .], the first and second order Frechet derivatives of ^ with respect to 
V can be obtained using the chain rule as 

^'(v)[v,]=^'{p)[p'{v)[v,]] (15) 

and 

^"{V)[VI,V2]= ^"{p)[p'{v)[v,lp'{v)[v2]] + ^'{p)[p"{v)[vuV2]l (16) 

where fi,f2 e S^. Equations (jlsj) and ((l6|) themselves can be easily expressed in 



terms of TZ{p) and its derivatives using (11) and (12). These equations will be 
used later as the key equations in finding the sensitivities in our parametric level 
set representation of p. 

3.2 Pixel Based Minimizations (Revisiting Traditional Level 
Set Methods) 

In the specific context of the shape-based inverse problems of interest here, the first 
and second order sensitivities of the cost function ^ with respect to the functions 
defining p(x) in ([T]) i.e., 0(x), pj(x) and Po(x), can be used to form a minimization 
process. Based on the order of the sensitivities available, first order optimization 
methods such as gradient descent or second order methods such as Newton or 
quasi-Newton techniques can be implemented. 

For simplicity in reviewing the current methods, we assume that Pi and Po 
are known a priori and only the shape (i.e., the zero level set of (p) is unknown 



see 25 , 28 , 76 for details on the recovery of both the shape as well as the contrast 



function). In an evolution approach it is desired to initialize a minimization process 
with some level set function 0o and evolve the function to find a cf) which minimizes 
^ . To take into account the concept of evolution, an artificial time is defined where 
the level set function at every time frame t > is rewritten as 0(x; t) and the zero 
level set of (/)(x; t) is denoted as dDt- A straightforward differentiation of 0(x; t) = 
with respect to t yields to the Hamilton- Jacobi type equation 

d(j) 



dt 



+ V(x;t) 







(17) 



for the points on dDt where ^(x; t) = dx/dt. To move the interface in the normal 
direction, \^(x; t) should be chosen as f (x; t)n(x; t) where w is a scalar speed func- 
tion and n = V0/|V0| is the unit outward vector on dDt. Incorporating this into 
the minimization of ^, the speed function for the points on dDt, denoted as v, 



can be chosen to be in the steepest descent direction of ^ which is 24 



v = -Me{{po-p^)n'{py[n{p)]}. 



(18) 



As (18) is only valid for x € dD, a velocity extension should be performed to ex- 



tend t; to f defined over the entire domain f2 and therefore capable of globally 
evolve the level set function 56 . Beside this classical level set approach, other 



ways of representing the speed functions and performing the minimization process 



are proposed 10,50 58 . For example, Hadj Miled and Miller 50 proposed a nor 



malized version of the classic speed function in the context of electrical resistance 
tomography. 

Some authors have also proposed using Newton type methods to update v at 
every iteration (e.g., see |To|[64||69]). Analogous to ^ and ([9]), in these problems 
the residual operator and the cost function are usually written directly as TZ^D) 
and ^{D), functions of the shape itself. Assuming the existence of the first and 
second order shape derivatives J5^, denoted as ^'{D)[ . ] and ^"{D)[., .], at every 
time step the Newton update v is obtained by solving [11] 

^"{D)[w,v] + ^'{D)[w] = Q yw^^D. (19) 

Here So is an appropriate Hilbert space such as L'^{dD), which may depend on 
the current shape (9,10 . Considering the general forms of the derivatives as (11) 
and (12), in a Gauss- Newton method the second derivatives of 71 are disregarded 
and (19) becomes 



^e{n'{D)[vln'{D)[w])s^+^e{n'{D)[wln{D))s^ = \fweSD. (20) 



Furthermore, to avoid ill conditioning, in a Levenberg-Marquardt approach (20) is 
regularized as 

^e(7^'(D)[i)],7^'(Z})H)§„+^e(7^'(Z})H,7^(Z}))s„+A(^},u;)§^ = ywe^D. 

(21) 



in which for A > 0, the equation is shown to be well-posed 10 . Similar to the 
previous approach, once v is obtained, a velocity extension is performed to result a 
globally defined v which can be used to update (p. However, although using second 
order methods can reduce the number of iterations in finding a minima, compared 
to gradient descent methods, they do not necessarily reduce the computation load. 
From an implementation perspective there are some concerns using the afore- 
mentioned methods. The gradient descent method usually requires many iterations 
to converge and performances become poor for low sensitivity problems l50j. Al- 
though using Newton and quasi-Newton methods to update the level set function 
increases the convergence rate, they are usually computationally challenging and 
for large problems and relatively finer grids, a large system of equations must 
be solved at every iteration. Also for both types of methods, there are usually 
added complications of the level set function re-initialization and speed function 
extension. The approach that we will put forth in the next section is capable of 
addressing these problems. It is low order and numerically speaking, the number 
of unknowns involved in the problem are usually much less than the number of 
grid points and hence allows us to easily use second order methods. Moreover, our 
proposed method does not require re-initialization of the level set function, speed 
function extension or even length-type regularization as a common regularization in 
many shape-based methods (e.g., see |24|[33|[50] ) and our level set function remains 
well behaved through the corresponding evolution process. 

4 A Parametric Level Set Approach 

As discussed earlier, in most current shape-based methods 0(x) is represented by 
function values on a dense discretization of x-space as part of a discretization of 
the underlying evolution equation or Newton-type algorithm. Consider now the 
level set function to be still a function of x but also a function of a parameter 
vector fi = (/ii,/i2, •■■^/^m) £ I^™- In this case we define the continuous Lipschitz 
function (p : Q x M"^ ^ M, as a parametric level set (PaLS) representation of D if 



for a c 6 M 

0(x, ^i) > c Vx e D 

0(x,^)=c y^edD (22) 

(/)(x, /x)<c Vxef2\D. 

In the PaLS approach we assume that the general form of 0(x, fi) is known and the 
specification of n can exphcitly define the level set function over the entire domain 
Q. In other words the evolution of required to solve the underlying inverse 
problems is performed via the evolution of /j,. An example of a PaLS function is 
a basis expansion with known basis functions and unknown weights and as will be 
shown later in this paper we considerably expand on this notion. We call fj, the 
PaLS parameters. 

To setup the problem using the PaLS approach, consider momentarily that 
Pj(x) and Po(x) are known a priori (later we will appropriately take away this 



restriction) . Based on ( 22 ) p is written as 



p(x, fx) = p,(x)if (0(x, ^) - c) + Po(x)(l - i7(0(x, /^) - c)). (23) 

Under this model, we now view ^^ in (|9| as a function of fi, i.e. ^(fJ-) '■ M™ -^ 
M. Therefore unlike the classic level set approach, the unknown is no longer the 
function 0, but a vector belonging to M'", where m is usually much smaller than 
the number of unknowns associated with a discretization of x-space. With this 
model and assuming (1) is sufficiently smooth with respect to the elements of 
H and (2) the discontinuous Heaviside function is replaced by a C^ approximation 
(e.g., [84]) denoted as Hrg, we now proceed to formulate a second order approach 



for the minimization of ^(fi). To begin, rewriting (23) with Hrg and taking a 
derivative with respect to yields 

^ = {p^-Po)Srgi<p-c), (24) 

where Srg{-) is accordingly the regularized version of the Dirac delta function (see 
examples in [73|[84]). Using the chain rule gives 

Op Op OS Oct) 



Now using (25) with (15) and (14), the gradient vector for ^ is 



dfXj 






dp 



^^e{n'{pnn{p)i{p,-po)6r,(<i>-c)^h,. 



(26) 
(27) 



We denote by J^(t^) the gradient of JF with respect to the parameter vector fj,. 
With this notation a gradient descent equation can be formed to evolve the PaLS 
function as 

(28) 



t^ 



(t+i) 



M 



(t) 






A'=M^ 



where A^*) > is the iteration step (m) and (28) is assumed to be initiahzed with 
some ^(°). Ahhough gradient decent is relatively simple to implement, it is known 
to be slow to converge and can suffer from difficulties associated with scaling of 
the parameters [21]. Moreover, the use of gradient decent fails to take advantage 
of one of the primary benefits of the PaLS idea; namely the ability to specify a 
level set function using a small (relative to a discretization of x-space) number of 
parameters. Under this model, it becomes feasible and indeed useful to consider 
higher order optimization methods such as Newton or quasi- Newton methods which 
are faster in convergence and robust with respect to sensitivity scahngs of different 
parameters p^. These methods usually use the information in the Hessian, which 
we now derive. To calculate the elements of the Hessian matrix for ^^ using (25) 
we have 



Q2p 



dfijdfik 



{Pi-Po){5rg{<p-c) 



dfijdfik 



+ (5' (0-c) 



d^j dfii 



-). 



(29) 



where S'^g(.) is the derivative of the regularized Dirac delta function. Based on (16) 
and (12) we have 






^"(pi^.^h^'ipi-"^" 



- dfij ' djik 



dnjdni 
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^^(^'(p)[|^],^'(p)[|^]k+^e(7^"(p)[|^,|^],7^(p))s. 
dfij dfik Ofij dfik 



+ ^e{n'(py[TZ{p)i 



Q2p 

dfijdfik 



h- (30) 



This equation is in fact the exact expression for the elements of the Hessian matrix 
for t^. However as mentioned earher, in methods such as the Gauss- Newton or 
Levenberg-Marquardt, to reduce the computation cost the Hessian is approximated 
in that the terms containing second order derivatives are disregarded. Following 



that approach here, (30) becomes 



^'^ -^e(7^'b)[|^],7^'(p)[|^])s„ (31) 



= Me{n'{p)[{p,-po)6rg{(l)-c)^ln'{p)[{p,-po)6rg{(l)-c)-^]h^. 

(32) 

We denote as H^(t^) the approximate Hessian matrix, the elements of which are 



obtained through (31). Having this in hand, a stable and faster converging PaLS 
function evolution can be proposed by solving the following Levenberg-Marquardt 
equation for ^(*+^) 

[h.(^)|,.,(. + AWi](m(*^^) - A^^*^) = -J.(^)|,.,(. t > 0. (33) 

Here A^*^ is a small positive number chosen at every iteration and I is the identity 



matrix. In fact referring to (31) we can see that the approximate Hessian matrix 
H;x(^) is a Gramian matrix and hence positive semidefinite. Therefore adding 
the small regularization term A(*)I would make the matrix at the left side always 
positive definite and hence the system of equations resulting the updates for ^ 
always has a unique solution. The left hand side matrix being strictly positive 
definite guarantees A/z^*) = /^(*+^) - ^i^^^ to be in the descent direction since we 

[AAi(*)]^[J^(^)] = -[A^W]T[H^(^) + AWI][Aax(*)] < (34) 

where the superscript ^ denotes the matrix transpose. More technical details about 
the implementation of the Levenberg-Marquardt algorithm such as the techniques 
of choosing A*^*) at each iteration based on trust region algorithms are available 
§|21 



m 



We now turn our attention to the determination of Pi(x) and Po(x). For sim- 
plicity here we assume that p(x) is piecewise constant as is often the case in work 



of this type 72 . This, in addition to the shape parameters, we need to deter- 
mine two constants defining the contrasts over the regions. We do note that the 
approach can easily be extended to consider other low order "texture models" as 
is done in e.g. |37]. As our primary interest in this paper is a new method for 



representing shape, we defer this work to the future. Under our piecewise constant 
contrast model, we denote the contrasts as pj(x) = pi and Po(x) = po, were pi and 



l-^ = Hr,ict>-c). (35) 



Po are unknown constant values. Following analogous equations as (26) and (30), 

these parameters can also be appended to the unknown PaLS parameters. The 

sensitivity of ^ with respect to these parameters can also be derived based on the 

fact that 

dp dp 

dpi dpo 

and the second order derivatives of p with respect to pi and po are zero. 

For the PaLS approach represented in this section, the intention is to keep the 
formulations general and emphasize the fact that this representation can formally 
reduce the dimensionality of the shaped based inversion. Clearly the expression 
d(j)ldfij depends on how the PaLS functions are related to their parameters. In the 
next section a specific PaLS representation is presented the efficiency of which will 
be later examined through some examples. 

5 PaLS Function Representation 

5.1 Adaptive Radial Basis Functions 

As pointed out earlier, appropriate choice of a PaLS function can significantly re- 
duce the dimensionality of an inverse problem. In this paper we are interested in a 
low order model for the level set function which will provide fiexibility in terms of 
its ability to represent shapes of varying degree of complexity as measured specifi- 
cally by the curvature of the boundary. This representation may allow for "coarse 
scale" elements capable of representing boundaries of low curvature with few ele- 
ments. Furthermore, it is desired to have finer grain elements capable of capturing 
higher curvature portions of the boundary such as sharp turns or corners. Such 
adaptability is desirable for a general PaLS representation since the characteristics 
must be present for well-posed inverse problems such as full view X-ray CT or 
even image segmentation where high fidelity reconstructions are possible. On the 
other hand, as we demonstrate in Section [6} for severely ill-posed problems, the 
availability of models with this parsimonious, but fiexible structure may allow for 
the recovery of geometric detail that otherwise would be not be obtainable form 
e.g., a traditional level set or pixel based approach. Writing a PaLS function as 
a weighted summation of some basis functions may be a reasonable choice here, 
where different terms in the summation may handle some desired properties about 



the shape. Here we focus specifically on the class of radial basis functions (RBF). 
We are motivated to concentrate on RBFs as they have shown to be very flexi- 
ble in representing functions of various detail levels. This flexibility makes them 
appropriate choice for topology optimization [79j , solving partial differential equa- 
tions 36 and multivariate interpolation of scattered data 32,81. Some examples 
of commonly used RBFs are Gaussian, multiquadric, polyharmonic splines and 
thin plate splines. More details about these functions and their properties are 
available in [SJ. 

Based on the statements made, consider the PaLS function 

mo 

0(x,a) = ^a,^(||x-x,||) (36) 

where ^ : IR+ ^ M is a sufficiently smooth radial basis function, a e M™" is the 
PaLS parameter vector and ||.|| denotes the Euclidean norm. The points Xj ^-^'s 
called the RBF centers. In an interpolation context, usually the centers are decided 
in advance and distributed more densely in regions with more data fluctuations. 
However, in a PaLS approach there may be limited information about the shape 
geometry and therefore more flexibility is required. Thus here we consider a more 
general PaLS function of the form 

mo 

<P{^,[(x,f3,x]) = Y.^MWA^-Xj)P), (37) 

i=i 

for which the vector of centers x = [Xi)X2r";Xmo] ^^^ ^^e dilation factors f3 = 
[/3i,/32,-")/3mo] are added to the set of PaLS parameters. Also in order to make 
the PaLS function globally differentiable with respect to the elements of /3 and x, 
similar to [I] a smooth approximation of the Euclidean norm is used as 



x||2 + t;2 VxeM", (3^ 



where t; ^^^ is a small real number. The use of (37) rather than (36) makes 
the PaLS function capable of following more details through scaling the RBFs or 
floating centers moving to regions where more details are required. To incorporate 
this into the optimization methods described, the sensitivities of (p with respect to 
the PaLS parameters are 

Hm^-XjW) (39) 



daj 



and 

Also considering Xj and a;^'^) to be the k^'^ components of Xj a-iid x as points in 
M", for k = l,2,---,n we have 

Clearly the sensitivities obtained are general and valid for any RBF in C^(]R+). In 
the next section, we consider a specific class of RBFs that we have found to be 
particularly well suited to the PaLS problem. 

5.2 The Choice of Compactly Supported Radial Basis Func- 
tions 

Usually the RBFs used in various applications involving function representations 
are C°° functions with global support [8] . However a different class of RBFs which 
have recently been under consideration are the compactly supported radial basis 
functions (CSRBFs) [sl]. These functions become exactly zero after a certain 
radius while still retaining various orders of smoothness. From a numerical point 
of view, compact support of the RBFs yields sparsity in the resulting matrices 
arising in the implementation of these methods and hence reduces the computation 
cost. This was recently the motivation to use these functions in simplifying level 



set methods 29 . Another interesting property of these functions is their local 
sensitivities to the underlying parameters ^I]. In other words, when a function is 
expressed as a weighted sum of CSRBFs, changing a term would not have a global 
effect on the function and only locally deforms it. 

Beside aforementioned advantages of the CSRBFs, our interest in this class of 
RBFs arises from their potential in reconstructing the shapes with a very small 



number of terms in a PaLS representation as (37). Furthermore, as will be ex- 
plained, this representation can involve shapes with corners and rather high cur- 
vature regions. For ■?/' > being a smooth CSRBF, lets denote every basis term in 



(l37j) as 

^,(x)=^(||/3,(x-x,)||^) (42) 



and call ipj a bump. Due to the compact support of these functions, for every two 
bumps tpj and ipk we can write 

supp{ipj + ipk) = supp(iljj) u supp(?/'fc). (43) 

For a real valued function i} defined over M" we define 

^c(^):={x:?9(x)>c}. (44) 



Clearly for c > 0, J^d"^) represents the interior of the c-level set of i). Based on (43 ) 
and the smoothness of the bumps, we obviously have that as c -» 0+, J^d'^Pj + V'fc) 
would tend to J^di^j) U-^(V'fc)- More generally for aj > 0, J^(Ej='i Q^j^j) tends to 
U^i-^(V'i) as (c/aj) -^ 0+. In other words, using relatively large CSRBF weights 
(as compared to c) can imply reconstruction of the shape through the union of a 
collection of floating balls of various radii. Moreover J^d'ipj - Oikipk) would tend to 
^c{,'>Pj)^^c{'>Pk) as c -* 0+ and ak -^ +oo. In a more general fashion, J^c{ciji'j-(yk4'k) 
tends to J^di^j)^ '^d'^k) as {c/aj) -^ 0+ and (ak/aj) -^ +cx). Therefore in this 
context, bumps with larger negative coefficients can yield holes or inflect the shape 
by excluding some portions of it. We would consider the two aforementioned 
properties of the CSRBFs as a "pseudo-logical" behavior of these functions. This 
property can result in rather high curvature geometries with a limited number of 
bumps. Besides high curvature regions, low curvature segments (e.g., an almost 
straight line in M? or a planar segment in M^) can be formed by interaction of two 
identical bumps but with opposite signs at their footprint intersections. Figure 
[ija sheds more light on aforementioned facts, and shows the interaction of two 
bumps, which for instance can represent a crescent with two rather sharp corners 
(a = -50), or a contour with a low curvature segment (a = -1). As a second 
example, a representation of a square using only 5 bumps is depicted in Figures 
HJb and [He. 



The most commonly used CSRBFs are those called Wendland's functions 81 
The smoothness and the compact support provided by Wendland's functions are 
of interest and hence we shall use them as the basis for our PaLS approach. Wend- 
land's functions follow the general form of 



vvM= n"'"" "^: («) 




when representing an RBF in M."^, with P„ ; being a univariate polynomial of degree 
[n/2j + 3/ + 1. In terms of smoothness, this class of RBFs belong to C^K A 






Figure 1: (a) Left: A close to zero (c = 0.01) level set of the function ipi + 01^2 for 
various values of a. The CSRBF used is the Wendland's function tpi^i (cf. Table fl]), 
the dilation factors are f3i = (32 = I and the centers are taken as xi = (^i; ^i) and 
;\;2 = (|, |). (b) Center: A PaLS function involved in representation of a square at 
a close to zero level set. Only 5 bumps are involved, (c) Right: The representation 
of the square 



Table 1: Compactly supported RBFs of minimal degree where 
(1 - r)+ = max{0, 1 - r} 



[n/2\ + 1 + 1 and 



Function 



V'n,i(r) = (l-r)fi((^+ 1)^ + 1) 
^n,2(r) = (1 - r)f 2((^2 + 4^ + 3)^2 + (3£ + g)^ + 3) 

'^^ 3(r) = (1 - r)f 3((£3 + 9£2 + 23£ + I5)r^+ 

(6£2 + 36£ + 45)r2 + (15£ + 45)r + 15) 



Smoothness 



^2 



derivation of these functions is provided in [8lj, from which we have listed the first 
few functions in Table [U 

In the next section we discuss the regularized heaviside function and explain 



how choosing an appropriate version of this function in (23), can pave the path for 
exploiting the pseudo-logical behavior of the bumps. 



5.3 Numerical Approximation of the Heaviside Function 



In a shape-based representation such as (23), solving the inverse problem numer- 



ically and making the evolution of the level set function possible requires using a 
smooth version of the heaviside function. A possible C°° regularization of H(.) 
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Figure 2: (a) Left: Two regularized versions of the heaviside function, (b) Right: 
Corresponding regularized delta functions 



denoted as Hi ^ is the one used in 116], as 



Hi^^(x) = -(l + — arctan( — )). 



(46) 



This function is commonly used in shape-based applications and specifically in the 
recent parametric representations of the level set function in image segmentation 
[5|[29]. Chan et al. in y^ have studied the characteristics of the resulting level set 
functions using Hi^^ as the evolution proceeds. In the context of the current shape- 
based problem for which Pi(x) = pi and Po(x) = po, and considering the zero level 
set, ([T]) reveals that in an evolution process (f) is likely to evolve towards a state that 
i/((/)(x)) = 1 for X e D and if ((^(x)) = for x 6 fi \ D. Referring to Figure [21a, one 
observes that when Hi,, is used as the regularized heaviside function, in order to 
have ifi^e(0) 2; 1 in D, the level set function should take rather large positive values 
(as compared to e) in this region. Analogously in fi \ D, is pushed to take rather 
large negative values. These constraints are implicitly imposed on the resulting 



level set function and specifically using CSRBFs in a PaLS approach as (37), the 



bumps are expected to distribute throughout Vt to form a level set function which 
takes rather large positive (or negative) values inside (or outside) D. 



An alternative choice of the regularized Heaviside function is the C^ function, 



1 x> e 

H2,e{x) = ] x<-e 

i + f + 2^sin(^) |a;|<e, 



(47) 



as proposed in 1841. It can be shown that if the nonzero level set c and e are 
appropriately set, using if2,e enables us to exploit the pseudo- logical behavior de- 
scribed previously. More specifically, for c > 0, and being a weighted sum of some 
bumps, and therefore compactly supported itself, we clearly want the points for 
which < to belong to VL\ D, i.e., to correspond to the intensity po- Using (23), 
this requires having H{(j)- c) = for < 0. Referring to Figure |2}a, this condition 
is satisfied when if2,e is chosen as the regularized heaviside function and -c < -e 
(or in general case \c\ > e as the required criteria). Practically, this use of if2,e takes 
away the implicit constraint imposed on the level set function using Hi,:, i.e., the 
bumps do not have to spread throughout f2 and may only concentrate inside and 
about the shape D to perform a higher resolution shaping. Figure |3] shows a typical 
shape representation resulted using Hi^ and H2^e highlighting the pseudo-logical 
property. 

Besides exploiting the pseudo-logical behavior, using i72,e can provide another 
advantage which reduces the dimensionality of the problem at every iteration, i.e., 
at every step of the evolution process, the cost function will be sensitive to a specific 
group of PaLS parameters. To further describe this behavior, consider the term 
5rg{4>-c)-r^ appearing in both the Jacobian and the approximate Hessian of ^ in 



dpi 



(27) and (32). Also for an evolution iteration, assume the superscript (t) indicating 



the value of every quantity at that iteration. If for a PaLS parameter such as ^ 



5rg{<P^'^-c) 



50(*) 



%o ^^0=40' 



Vxefi, 



(4J 



then using either one of the minimization schemes as (28) or (33) yields 



M 



,(*+!) _ ,,(*) 



Jo 



M 



Jo 



(49) 



The reason for this result is clear for the gradient descent scheme (28), based 



„(*) 



on the fact that using (48) in (27) causes the j^^ element of gradient vector to 



vanish and hence /r- remaining unchanged at the corresponding iteration. For 



Jo 



the Levenberg-Marquardt scheme (33), using (48) in (32) causes all the elements 
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Figure 3: (a) Left: A typical PaLS function resulted using Hi,, with 65 bumps and 
considering the zero level set. (b) Right: A typical PaLS function resulted using 
if2,e with 30 bumps and considering the c level set 



in the j*^ row of the approximate Hessian matrix to vanish, which results the 
corresponding equation 

A(/^f)-/^.;^)=0, (50) 



again equivalent to (49). Therefore, for either one of the minimization schemes, if 



rt48| holds, the PaLS parameter ^ will stay unchanged in that iteration. We here 
describe a common case that (48) holds during the evolution process: 

By using if2,e as the regularized level set function, the corresponding regularized 
delta function 52^e will be compactly supported (as shown in Figure |2]b), hence 
'52,e(0 - c) is only nonzero for c-e<0<c + e. On the other hand, based on the 
PaLS approach presented in this paper using CSRBFs, for /ij being any of the 
PaLS parameters a,-, /?,• or x) corresponding to the bump ■?/',•, we have 



supp(^) £supp(?/^j). 



(51) 



This fact is easily observable in (39), (40) and (41), where the related derivatives 



can only have nonzero values in snpp{ipj). Therefore, if at some iteration and for 
a bump ipj , 

supp (52,e(0 - c)) n supp(V^jJ = 0, (52) 

then in that iteration we have 

S2M-c)^-0 (53) 

and therefore the PaLS parameters corresponding to ipj will stay unchanged in 
that iteration. Figure |4] illustrates this phenomenon, showing a PaLS function 
composed of 6 bumps at some iterations. For 5 of the bumps used, the corre- 
sponding support does intersect the region supp(52,e(0 - c)), and therefore their 
corresponding parameters have the potential to change at this state of the PaLS 
function. However, a bump denoted as ijjj^ , does not intersect supp ((52,e(0-c)), and 
the underlying PaLS parameters do not need to be considered in that iteration. 
This approach is similar to the narrow-banding approach in traditional level-set 
methods (2,59,66 , where the values of the level set function are only updated on 



a narrow band around the zero level set and hence reducing the computation load. 
In our approach, however, this band is the points for which c-e < (p < c + e and the 
bumps which do not intersect with this band do not evolve at the corresponding 
iteration and hence their corresponding parameters are not updated. 

In the next section, through a number of examples drawn from a wide range 
of applications, we will show the superior performance of the proposed method 
specifically exploiting the pseudo-logical behavior of the CSRBFs. 

6 Examples 

In this section, we examine our proposed method for three different inverse prob- 
lems, namely electrical resistance tomography. X-ray computed tomography and 
diffuse optical tomography. The examples are simulated for 2D imaging and the 
results are provided in each section. Throughout all the examples, Q denotes the 
region to be imaged and D denotes the shape domain as stated in the previous 
sections. 
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Figure 4: (a) Left: A typical PaLS function composed of 6 bumps, and the c ± e 
level sets. (b) Right: The bumps with active evolving parameters 

and the frozen bump 



6.1 Electrical Resistance Tomography 

As the first example we consider electrical resistance tomography (ERT) catego- 
rized as a severely ill-posed problem 22,31 . The objective of this problem is the 



reconstruction of the electrical conductivity within a region of interest based on the 
potential or current measurements performed at the periphery of that region. Such 



reconstructions may be applicable in various areas such as medical imaging 17 



geophysics 23 and environmental monitoring 19 



For many geophysical applications the underlying physical model describing 
the DC potential m(x) inside Q in terms of the conductivity cr(x) and the current 
source distribution s(x) is 



V • (o-(x)Vm(x)) = s(x) in i7 
du 



(54) 



ov 

M = 



on dVtn (= dVt, 

on dVt^ = dVt \ dVtr^ 



where u denotes the outward unit normal on dVt and dVtn and dVL^ correspond 
to Neumann and Dirichlet boundaries. In many of the applications (e.g., see 



45,54,65,75]) the Dirichlet boundary condition is imposed as an approximation 



to the potential in regions far from the actual imaging region, and is used here for 



simplicity. For arbitrary distributions of the conductivity, (54) is usually solved nu- 

However, 



merically by means of finite element or finite difference methods |69j , j83j 

as the main focus of the paper, we concentrate on piecewise constant conductivity 

distribution as cr(x) = cXi for x e D and o"(x) = cTq for x € f2 \ D. 

For the inverse problem, the sensitivities of the measurements to perturbations 
of the conductivity (in our approach the perturbations of the PaLS parameters) 
are required. For s(x) - 5(x-Xs), i.e., a point source current at x^ e fi, we denote 
by Ms(x) the resulting potential over the domain fl and consider the measured 
potential at x^^ e i7 as 

Uds= / Ms(x)(5(x-Xrf)dx. (55) 

Jn 

The variation of Uds resulting from a perturbation 6a in the conductivity (i.e., the 
Frechet derivative of the measurements with respect to the conductivity) can be 



then expressed as 62 , 68 



duds 
da 



[6a] = f 
Jn 



6a VMs • VWd dx. 



(56) 



where Ud is the adjoint field that results from placing the current point source 
at X(i. To express the inverse problem in a PaLS framework, we consider M.(.) as 
the nonlinear forward model mapping the conductivity distribution into a vector of 
voltage measurements u obtained by performing M experiments, having a different 
point source position at each experiment and making A^^ potential measurements 
for £ = 1,2,---,M. Having the residual operator 71(a) = M.{a) - u and using (56), 



the Frechet derivative denoted as 'R''{a)[.] can be considered as a vector consisting 
of M sub-vectors Tl'^(a)[.], structured as 



n',{a)[6a] 



J^ 6a Vue ■ VMi dx 
In ^^ VMf • Vu%^ dx 



(57) 



Here ui denotes the potential in the i^'^ experiment and u^ denotes the adjoint 
field corresponding to the i^'^ experiment resulted from placing the current point 
source at the i^^ measurement point. Having Tl'(a)[.] in hand, one can obtain the 



PaLS evolution through using (26) and (32) in (33) 



For the purpose of this example, we model the electric potential within the 
box Q = [-3,3] X [-3,0] all dimensions in units of meters in x - y plane. Here 
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Figure 5: (a) Left: The 2D modelling region for ERT. The darker interior region 
is the imaging region surrounded by the sources and detectors. The dashed lines 
correspond to dipole nodes used in every experiment (b) Right: The gray region 
shows the shape to be reconstructed in the ERT problem. The dots with "+" and 
"-" signs correspond to the centers of positive and negative weighted bumps in the 
initial state of the problem. The black contour is the resulting c-level set of the 
initial PaLS function. 



dfln corresponds to the top surface {y = 0) and dfld corresponds to the sides and 
bottom boundaries as shown in Figure [5ja. The region to be imaged is the square 
Q' = [-0.5, 0.5] X [-1, 0], where 30 points are allocated to the sensors placed equally 
spaced on the top and the sides of this region (shown as small circles in the figure) . 
A total of M = 40 experiments are performed and in every experiment two of the 
sensors are used as a current dipole (superposition of two point sources with oppo- 
site signs) and the remaining 28 sensors measure the potential at the corresponding 
positions. The dipole sources are chosen in a cross-medium configuration, where 
the electrodes corresponding to every experiment are connected with a dashed line 
as again shown in Figure [5ja. With this configuration we try to enforce electric 
current flow across the medium anomalies and obtain data sets more sensitive 
to shape characteristics. For the simulation, we use the finite difference method 
where Q is discretized to 125 grid points in the x direction and 100 points in the 
y direction. The gridding is performed in a way that we have a uniform 75 x 75 



grid within Q' (excluding the boundaries containing the sensors) and the exterior 
grids hnearly get coarser as they get further from Q' in the x and y direction. The 
forward modelhng is performed over the whole collection of grids in Q, while the 
inversion only involves the pixels within fi' (known as active cells |60]). 

The shape to be reconstructed is shown as the gray region in Figure |5}b. This 
shape is a threshholded version of a real scenario and is of particular interest here 
because it has a concavity facing the bottom where there are no measurements 
performed. The values for the anomaly conductivities are ai - 0.05 Sm" and 
(To - 0.01 Sm~^ and the conductivity value used for ^2 \ f2' is the same as the true 
value of o"o in all our inversions. In the inversions we consider the data u to be 
the measured potentials generated by the true anomaly with 1% additive Gaussian 
noise. For the shape representation we use if2,e with e = 0.1 and we consider the 
c = 0.15 level set. 

For the PaLS representation (36) we have mo = 40 terms and the Wendland's 
function ipi^i is used as the corresponding bump. For the initial PaLS parameters, 
we consider a random initial distribution of the centers x , within the square 
[-0.4,0.4] X [-0.8,0]. The weighting coefficients are initialized as aj = ±0.2 where 
the centers of alternatively positive and negative weighted bumps are shown with 
" + " and " - " signs in Figure [5jb. The positive initial values of aj are taken 
slightly bigger than c to have some initial c-level sets. The purpose behind having 
alternative bump signs is to have the narrow-band supp((52,e(0 ^ c)) cover various 
regions of Q' and increase the chance of initially involving more bumps in the 
shaping as explained in previous section and Figure |4j The dilation factors are 
taken uniformly to be Pj = 4, as an initialization to make the support radii of 
the bumps small enough for capturing details and more or less large enough to 
carpet the region Q'. Our intention for this initialization of the PaLS parameters 
is to provide a rather simple, reproducible and general initialization. The stopping 
criteria in the reconstructions is when the norm of the residual operator reaches 
the noise norm (known in the regularization literature as the discrepancy principle 

Figure [6ja shows the shape reconstruction using the PaLS approach, and as- 
suming the values ctj and cTq are a priori known. Figure [6jb shows the result of the 
same problem using a typical Gauss-Newton approach with the level set function 
defined as a signed distance function over the pixels and a smoothing regulariza- 
tion added as explained in l69j. This algorithm only considers shape reconstruction 
(i.e., (Tj and (Tq are considered known) and it is initialized with the same contour 
as the initial c-level set contour of the PaLS approach. As the results in Figure 
|6|a and[6jb show, the PaLS approach performs well in reconstructing major shape 
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Figure 6: (a) Left: The result of reconstructing the shape using PaLS approach 
after 26 iterations, the centers Xj ^iid their corresponding weight signs are shown, 
(b) Center: Result of using the traditional level set method after 39 iterations (c) 
Right: The result of reconstructing both the shape and the binary conductivity 
values after 32 iterations 



characteristics, while the traditional level set approach fails to provide a good re- 
construction in low sensitivity regions close to the bottom and does not capture 
the concavity. Figure [6]c shows the result of reconstructing both the shape and the 
anomaly values using the PaLS approach, where this time the PaLS evolution takes 
slightly more iterations (32 iterations verses 26), but the resulting reconstruction 
still well represents the shape. The initial values used for the conductivity val- 
ues are a\ = 0.01 Sm~^ and ao = 0.005 Sm~^ and the final resulting values are 

^32^ —1 f32^ —1 

(Tj - 0.056 Sm and ao = 0.010 Sm , which show a good match with the 
real values. To illustrate the behavior of the PaLS function, in Figure [7] we have 
shown the initial and final PaLS functions for the shape only reconstruction. Also 
to compare the convergence behaviors of the PaLS approach and the traditional 
level set approach, in Figure [8] we show the residual error through the evolution 
steps for both methods. Using the PaLS approach the stopping criteria is met 
after 26 iterations while traditional level set method reaches a local minima after 
39 iterations (the updates after 39 iterations become so small that it stops evolving 
further) . 





Figure 7: (a) Left: The initial PaLS function (b) Right: The final PaLS function 
for the shape-only reconstruction of Figure 16} a 



6.2 X-ray Computed Tomography 

As the second example and a mildly ill-posed problem, we consider X-ray Com- 
puted Tomography (CT) (20]. In this imaging technique, X-ray photons are trans- 
mitted through the object of interest and the intensity of transmitted ray is mea- 
sured at the boundaries to reconstruct the object's attenuation coefficient. The 
contrast between the attenuation characteristics of different materials can provide 
structural information about the object being imaged. X-ray CT is among the 



most well known methods for imaging the body tissue in medical applications 35 



For an X-ray beam transmitting along a line Ck in the tissue, the photon 
intensity Xk measured at the detector side of the line can be written as 



Xk= I Ik{£)exp[- / a(x,£')dx) dS 



(5^ 



where a(x, £^) denotes the attenuation coefficient, in general as a function of the 
position X and the energy of the incident ray £, and Ik(£) denotes the incident 
ray energy spectrum. In case of a monoenergetic beam as Ik{£) = Io,kS(,S ~ ^o), a 
measured quantity related to the photon intensity may be defined as 



Uk 



logCy-^) = / a(x)dx. 

-'0,fc -'•Cfc 



(59) 
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Figure 8: Residual error reduction through the iterative process for the shape only 
reconstruction, using PaLS approach and the traditional level set method 



Equation (59) simply relates the measurements to the Radon transform of the 



attenuation coefficient a(x) in a monoenergetic scenario. The quantities Uk are 
actually what is considered as the data in CT imaging. 

The Frechet derivative of the CT measurements with respect to the attenuation 
coefficient is expressed as 

^ [6a] = f 6a dx. (60) 

da JCk 

Considering u to be the set of CT data collected along different paths £k, for 

k = 1,2,---A^, and A4{a) as the forward model mapping the attenuation to the CT 

data set, based on (60) the sensitivity of the residual operator 71(a) = M.{a) - u 



with respect to a perturbation in a can be written as 

J^ 6a dx 



n'(a)[6a] 



[r 6a dx 



(61) 



which we need for the PALS evolution process. 

We consider 2D imaging over a square of 2 m x 2 m, i.e., Vt = [-1, 1] x [-1, 1] 
in the x - y plane, as shown in Figure lOla. The region outside Q is assumed to 
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Figure 9: (a) Left: The 2D X-ray CT imaging setup (b) Right: The gray region 
shows the attenuation shape to be reconstructed. The dots with "+" and "-" signs 
correspond to the centers of positive and negative weighted bumps in the initial 
state of the problem. The black contour is the resulting c-level set of the initial 
PaLS function 



have zero attenuation. The X-ray beams are emitted as parallel beams and the 
measurements are performed through an equally spaced linear array of 34 sensors, 
the vertical axis of which makes an angle Q with the x axis. For a full view X-ray CT 
we have 6o <6 <9o + tt, while in a limited view scenario, and hence a more ill-posed 
problem, 6 varies in a smaller angular domain. The shape to be reconstructed is 



shown in Figure 9 b, which is brought here from 72 as a rather complex shape to 



examine the flexibility of the PaLS approach. The region fl is discretized into 64x64 
uniform grid. The binary attenuation values are a^ = 2.5 cm^^ and Oo = 1 cm^^. 
For the purpose of imaging, the forward model measurements are performed at 
every 1 degree angle. For the PaLS representation, we use mo - 50 bumps with the 
centers Xjy distributed randomly as shown in Figure [91 b. For this example we use 
slightly more bumps due to the better posed nature of the problem (at least in the 
full view case) and the more complex shape. Again to roughly carpet the domain 
Q with the bumps we uniformly take the dilation factors to be Pj = 2.5. The other 
PaLS settings used are similar to those of the ERT example in the previous section. 
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Figure 10: (a) Left: The result of reconstructing both the shape and the anomaly 
values in a 5% noise case using the PaLS approach, convergence achieved after 
20 iterations (b) Right: Result of using the same data set with traditional level 
set method to recover shape only, assuming attenuation coefficients are known, 
convergence achieved after 14 iterations 



Figure 10 a shows the result of reconstructing both the shape and the atten- 



uation values for a full view experiment where < ^ < vr. The synthetic data are 
obtained from the true attenuation map and by adding 5% Gaussian noise to the 
measurements. The initial values used for the attenuations are cr^ - 1.5 cm"^ and 
do - 0.5 cm~^ and the discrepancy principle stopping criteria is met after 20 iter- 
ations resulting the corresponding shape and the final values oi oq = 2.494 cm^^ 
and OTo = 0.997 cm^^, which are very well matched with the real quantities. We 
should mention here that when the attenuation values are assumed to be known, 
the number of iterations for shape reconstruction is only 12. The same data set is 
used to reconstruct the shape only, using the traditional level set method used in 
the ERT example and initialized with the same contour as the c-level set contour 



shown in Figure |9jb. The result of this reconstruction is shown in Figure 10 b. 
Due to the better posed nature of this problem, the pixel based level set method 
also provides a successful reconstruction for the geometry only. Note, however that 
the PaLS approach is solving a more challenging problem (i.e., reconstructing the 
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Figure 11: (a) Left: The result of PaLS shape reconstruction in a 1% noise case 
to show the flexibihty of following shape details due to the pseudo-logical prop- 
erty, convergence achieved after 42 iterations (b) Right: The corresponding PaLS 
function 



attenuation values as well as the shape), thanks to the pseudo-logical behavior of 
the bumps, the resulting shape follows the same detail levels as a pixel based level 
set function. This shaping capability of the PaLS is more highlighted when we use 



a less noisy data (1% Gaussian noise). The result is shown in Figure 11 where the 
convergence is achieved after 42 iterations, and the resulting contours follows the 
true shape in high level of details. Finally as a more challenging and very ill-posed 
problem, we now consider a limited view scenario where the angular domain is lim- 
ited to 7r/4 < 9 < 37r/4, and fewer rays cross the anomalies. Two percent additive 
Gaussian noise is added to the synthetic data obtained from the true attenuation 



map. With the same problem setting as the full view case. Figure [12| a shows the 
result of reconstructing the shape using the PaLS approach. However, as shown 



in Figure 12 b, the traditional level set method applied to this problem given the 



attenuation values, fails to provide a complete reconstruction and stops further 
shape enhancement after reaching a local minima. 
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Figure 12: (a) Left: The PaLS shape reconstruction in a hmited view X-ray CT, 
convergence achieved after 49 iterations (b) Right: The performance of the tra- 
ditional level set method for the same data set, reaching a local minima after 32 
iterations 



6.3 Diffuse Optical Tomography 

As the third application we consider diffuse optical tomography (DOT). In this 
imaging method data are obtained by transmitting near-infrared light into a highly 
absorbing and scattering medium and then recording the transmitted light. As the 
inverse problem of interest the photon flux is measured on the surface to recover 
the optical properties of the medium such as absorption and/or reduced scattering 
functions. A well known application of this method is breast tissue imaging, where 
the differences in the absorption and scattering may indicate the presence of a 
tumor or other anomaly [T] . For given absorption and scattering functions a number 
of mathematical models have been proposed in the literature to determine the 
synthetic data (photon flux) Is] . We focus on the frequency-domain diffusion model 
in which the data is a non-linear function of the absorption and scattering functions. 
Consider Q, a square with a limited number of sources on the top and a limited 
number of detectors on either the top or bottom or both. We use the diffusion 
model in [s] modified for the 2D case where the photon flux ^(x;^) is related to 



input s{'x;u) through 



u 



- V ■ /3(x) Vm(x; lu) + a;(x)M(x; u) + i—u{yi] uj) = s(x; w), (62) 



V 



with the Robin boundary conditions 

M(x;a;) + 2/3(x)^^^|^ = 0, ^edQ. (63) 

Here, /3(x) denotes the diffusion, which is related to the reduced scattering func- 
tion //^(x), by /3(x) = l/(3/i^x)), a;(x) denotes absorption and d/du denotes the 

normal derivative. We also have i = v-T, to as the frequency modulation of light 
and V being the speed of light in the medium. Knowing the source and the func- 
tions q;(x) and /3(x), we can compute the corresponding u(x;ci;) everywhere, in 
particular, at the detectors. 

As the inverse problem, we consider a case that the reduced scattering function 
is known and we want to reconstruct the absorption a(x) from the data. Again 
for this problem we consider a point source s(x; u) = 5(x - x^) for x^ e dQ, which 
results in a photo flux Ms(x; u) in Q. For a measurement at Xf^ € i7 as 

Uds= / Ms(x;u;)(5(x-Xrf)dx, (64) 

Jn 

the variations with respect to the variations in the absorption can be written as [s] 

duds 



da 



■[6a] = / (5q; Ms(x;a;) M(i(x;cj)dx, (65) 

Jn 



with Ud being the adjoint field caused by having the point source at x^^. Identical to 
the notation used for the ERT problem, consider 71(a) = M.(a) - u as the residual 
operator where the data vector u contains the measurements of A^^ detectors from 



M experiments. Based on (65) the i^"- block oiTl'(a)[.] corresponding to the i^'^ 
experiment is 

(J^ 6a ue u{ dx 
; I , (66) 

with U£ denoting the photon flux in the i^'^ experiment and u^ being the adjoint 
field. For the purpose of this example we consider a very ill-posed problem where 
Q corresponds to a 5cm x 5cm imaging region, and there are only 8 point sources 
on the top and 8 detectors at the bottom. In every experiment only one of the 
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Figure 13: The 2D DOT imaging setup with the sources and detectors setting. 
The gray region shows the absorption shape, while 2% Gaussian noise is added to 
the absorption as a heterogeneity. The dots with "+" and "-" signs correspond 
to the centers of positive and negative weighted bumps in the initial state of the 
problem. The black contour is the resulting c-level set of the initial PaLS function 



sources is on and the measurements are performed at the bottom detectors. We 
have /is(x) = 6 cm~^ throughout Q. The measurements are performed at DC mode, 
25 MHz and 50 MHz. The absorption binary distribution is shown in Figure 13 
with values ctj = 0.015 cm^^ and ao = 0.005 cm^^. To make the problem more 
challenging we made the absorption distribution slightly heterogeneous by adding 
2% white Gaussian noise to it. For the PaLS setting, due to the very ill-posed 
nature of the problem we use relatively fewer bumps, i.e., mo = 20. Similar to 
the previous two examples the centers Xj corresponding to positive and negative 
weighted bumps are taken randomly inside Q as shown in Figure 13 Again to 



have reasonable initial support radius for the bumps the dilation factors are set 
to be Pj = 80. Other PaLS settings are similar to the previous examples. The 
forward model is solved using finite difference method by discretizing fi to 50 x 50 



grid points. Figure [T4| a shows the result of our shape-only reconstruction after 152 
iterations, when the true absorption map is used to generate the data and 0.1% 
Gaussian noise is added to it. In case of not having any noise in the data (but still 
considering the heterogeneity in the absorption), the results after 200 iterations are 





Figure 14: (a) Left: The PaLS shape-only reconstruction, stopping criteria 
achieved after 152 iterations in case 0.1% data noise and 2% heterogeneity (b) 
Right: The resuhs of shape reconstruction in case of only heterogeneity in the 
absorption. The evolution is manually stopped after 200 iterations, showing the 
results 



shown in figure [14} b. Although the current problem settings make it very ill-posed 
and extremely hard for shape reconstruction, the PaLS approach shows reasonable 
reconstructions, providing important details. For this very ill-posed problem, using 
the traditional level set method we failed to generate useful reconstructions. 



7 Concluding Remarks 



In this paper we proposed a parametric level set approach to shape based inverse 
problems. The basic formulation of the problem is kept general, emphasizing the 
fact that numerically, an appropriate parametrization of the level set function is 
capable of reducing the problem dimensionality and intrinsically regularizes the 
problem. Through such modelling the level set function is more likely to remain 
well behaved, and therefore there is no need to re-initialize as in traditional level 
set methods. Moreover, based on the fact that the number of underlying param- 
eters in a parametric approach is usually much less than the number of pixels 
(voxels) resulting from the discretization of the level set function, we make a use of 



a Newton type method to solve the underlying optimization problem. We specif- 
ically considered using compactly supported radial basis functions, which used as 
proposed provides two main advantages besides being parametric; first, the pseudo- 
logical property allows for recovery of wide range of shapes; and second, implicit 
narrow-banding for the evolution. 

Although in this paper, we only considered compactly supported radial basis 
functions as the bump functions with circular support, a combination of sufficiently 
smooth bumps with various types of supports may be considered in applications 
where there is more prior information about the shape and its general geometric 
structure. An efficient classification of such functions for various types of shapes 
and problems (i.e., forming a basis dictionary) would be a future direction of re- 
search. Although in the examples presented in this work the number of RBFs where 
coarsely chosen, we are currently developing some ideas to adoptively determine 
the number of basis elements using a sparse signal processing approach. Clearly 
this matter is beyond the scope of this paper and is planned to be presented as 
a future work. Moreover, for sake of simplicity in this paper we only considered 
2D problems, whereas the efficiency of a parametric approach is more pronounced 
for 3D shape reconstructions where the contrast between the number of voxels 
and the number of PaLS parameters in a parametric approach is more significant. 
Efficiently modelling 3D scenarios via both parametric shape representation and 
appropriate texture and heterogeneity models is an important future direction and 
a continuation of the current work. 
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